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0\ : Abstract 

The origin of self-organized criticality in a model without conservation 
. law (Olami, Feder, and Christensen, Phys. Rev. Lett. 68, 1244 (1992)) is 

QQ ' studied. The homogeneous system with periodic boundary condition is found 

o' 

t^^ ' to be periodic and neutrally stable. A change to open boundaries results in 



o 



the invasion of the interior by a "self-organized" region. The mechanism for 



C^ ' the self-organization is closely related to the synchronization or phase-locking 

a" 

of the individual elements with each other. A simplified model of marginal 



oscillator locking on a directed lattice is used to explain many of the features 
in the non-conserved model: in particular, the dependence of the avalanche- 



^ ' distribution exponent on the conservation parameter a is examined. 
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The phenomenon of the self-organized criticahty (SOC) |]T| is characterized by sponta- 
neous and dynamical generation of scale-invariance in an extended non-equilibrium system. 
One of the key issues in this field is to identify the mechanism(s) of SOC. It has been con- 
jectured that conservation laws or special symmetries are necessary 0. Conservation laws 
certainly are of great importance in "sandpile" models |I|J^J^, where the scale invariance 
can be shown to follow from a local conservation law (sand grains are conserved except 
at the boundaries of the pile) P]. In this sense, the origins of long range correlations in 
SOC systems with conservation are well understood, though not all exponents have been 
calculated analytically. However, models have been constructed ^^ that have no apparent 
conservation law, and yet display a power-law distribution of avalanche sizes. Of particu- 
lar interest is the model proposed by Olami, Feder, and Christensen (OFC) fl-pl which, 



they argue, models earthquake dynamics. The OFC model is very similar to the conserved 
sandpile model [|[] , but it has a parameter which defines the degree of conservation. In this 
paper, we study in detail the OFC model and we find that the self-organization is due to 
synchronization or "phase-locking" - a mechanism very different from that in the conserved 
models. 

In the OFC model, dynamical "height" variables hi are defined on sites i of a square 
lattice. The hi increase at unit rate until h = 1 a.t some location. The site j where hj = 1 
is considered to be unstable and will "topple". The rule of toppling is that when hj > 1, 
then hj — > and h^ ^ hj. + ahj, for all k neighboring j. The toppling on site j may 
cause its neighbors to become unstable {h^ > 1) and to topple. This procedure is repeated 
until all sites are stable {hi < 1 everywhere). The magnitude of the avalanche is given by 
the total "energy" dissipated in the process, i.e., the total change in J2ihi. The avalanche. 



which happens instantly on the time scale of driving |T3|, is then followed by growth. The 
parameter a is the measure of conservation (of /I's). When a = 1/4, the model is conserved 
and it is in the same universality class as for the BTW model |I|,Pjr^. We study here the 
non-conservative case < a < 1/4. 

We note that there is an ambiguity in this model. After some number of avalanches, it 



will occur that more than one site will have exactly the same height, as more than one site 
may topple to a height of exactly zero in a single avalanche. This can lead to two neighboring 
sites toppling simultaneously; the above procedure does not define the result of such events. 
We have modified the rules in several ways, e.g. by adding a very small amount of random 
noise to ensure that no two sites have the same height, and they all give the similar results. 

As most avalanches are small, and do not change the hi at many sites, we do not need to 
scan the whole lattice after each avalanche in order to determine the most unstable site. We 
use a tree structure to keep track of the highest values of hi in order to determine avalanche 
trigger sites. Using this technique, we have simulated up to more than 10^'^ avalanches for 
single systems. 

We find that the system with periodic boundary conditions quickly reaches an exactly 
periodic state |T^, with a unique period in the slow time. It has been noted [§,|l3] that in 



the case of periodic boundary conditions, the avalanche size distribution function drops very 
quickly with size. In the case of our modified model, which prevents two sites from having 
the same value of h, all avalanches consist of the toppling of exactly one lattice site, after 
a brief transient time. In a periodic state, the /ij's take turns to topple, one by one. The 
height h decreases by one each time a site topples and increases by 4a due to the toppling 
of the four neighbors when they reach h = 1. Thus the period of all these periodic states 
is 1 — 4a in the slow time variable, so that the slow "growth" is balanced exactly by the 
dissipation due to toppling. These periodic states are highly degenerate and neutrally stable 
in the sense that a typical small perturbation of the height at a single site in a periodic state 
is still a periodic state. They are similar to the neutrally stable periodic states in coupled 
oscillators [1^. There is a continuous set of periodic states in the attractor, with measure 



(1 — 4a) ^ in the initial phase space, where V is the system volume. 

Any inhomogeneity, such as a change in boundary conditions, destroys such simple peri- 
odic states. When the boundary conditions are open, the system can no longer have period 
1 — 4a, as the boundary sites have 3 neighbors (we study a system that is open on one axis, 
with the other directions periodic). Initially, the interior sites quickly converge to a nearly 



periodic state and topple with period 1 — 4a, but the boundaries are aperiodic. At longer 
times, the aperiodic region invades the periodic interior, as shown in Fig. [1|. This invasion, 
which destroys the periodicity in the interior and builds up long-range correlations, occurs 
by a mechanism similar to oscillator locking, as we describe below. The interface between 
the two regions is well defined on scales larger than one or two lattice constants. The inva- 
sion distance appears to have a power law dependence on time, y(t) ~ t'^*-"-' as shown in Fig. 
I, with (3 = 0.23 ± 0.08,0.63 ± 0.08 for a = 0.07,0.15, respectively. We see such invasion 
occurring even for values of a < 0.05, though f3 appears to be quite small, so that the 
invasion is extremely slow. This suggests that the transition to non-SOC behavior claimed 
by Olami, Feder, and Christensen |^ may only be apparent, due to the finite time of the 
simulations; we note that the time for complete invasion of a 128^ system with a = 0.07 is 
greater than 10^° avalanches. In the limit of long times, when the invasion crosses the whole 
sample, the distribution P{s; a) of avalanches of size s is a power-law, with 

P(s;a)~s-"("). (1) 

Consistent with OFC, we find r(a) = 3.2 ± 0.1,2.3 ± 0.1 for a = 0.07,0.15, respectively. 
The avalanches are not uniformly distributed in the system; the typical avalanche size grows 
with distance from the edge. 

The infiuence of inhomogeneity, which leads to a breaking of the 1 — 4a; periodic state, is 
quite strong. Even a single defect, where h is set to zero for all time, destroys the periodic 
state, and leads to a power law distribution of avalanche sizes. 

The spatial distribution of length scales is apparent in Fig. |l|(d), which shows an example 
of a configuration between avalanches in the steady state. Near the boundaries, the hi 
have only short range correlations, but this correlation length grows with distance from the 
boundary. The typical avalanches in the interior, though, still involve the toppling of only 
one lattice site. Infrequently, an avalanche is triggered near the boundaries, which penetrates 
into a longer-range correlated region, giving a large avalanche. We define a toppling rate 
r{y) as a function of distance y from the boundary, which gives the inverse of the mean time 



between topplings. At the boundaries, where the sites have only three neighbors, r{y) is 
smaller than the interior, where r{y) -^ (1 —Aa)^^ as y ^ oo. The toppling rate differential, 
defined as 6r{y) = (1 — Aa)^^ — r{y) is found to behave as 

6riy) ~ y-^ (2) 

with 1] = 3.2 ± 0.6, 1.8 ± 0.2 for a = 0.07, 0.15, respectively. Let ht{y) be the average height 
just before toppling, R{y) = (1 — Aa)r{y)ht{y) is then the dissipation rate and -R(y) ^ 1 as 
y ^ oo. It can be shown that in the steady state SR{y) = 1 — R{y) behaves as SR{y) = 



exp{—yJ ). Thus the dissipation rate is rather uniform except within a boundary layer 



of thickness J i^^^ - The power-law behavior of Eq. @j must be compensated by a power-law 
in 5ht{y) = ht{y) - 1: 

6htiy) ~ y-". (3) 

Note that the toppling height can be larger than 1 only when the toppling is triggered by 
neighboring sites and hence is part of an avalanche which involves more than one site. The 
picture we get from Eqs. (H) and (^ is that although the toppling happens less frequently 
as we move towards the boundary a larger portion of it is due to multi-site avalanches. 

In order to gain some insights on the build-up of long range correlations in the inho- 
mogeneous system, we consider a system which has only two sites: hi and /i2. Let us first 
consider the homogeneous case: both hi and /i2 are driven with unit rate and when one 
of them reaches the value one it topples. The rule of toppling is that: if /ii(2) > 1, then 
^2(1) -^ ^2(1) + tt^i(2) and /ii(2) -^ 0. This simple system is completely integrable. It has a 
continuous set of periodic states which are marginally stable. To illustrate its dynamics, we 
construct a Poincare map. Denote hi{n) to be the value of hi right after the nth toppling 



of /i2. It is easy to show that |T^ 



hi{n + l) = < 



hi{n) a < hi{n) < 1 

l + a-ahi{n) 1 <hi{n) <l/a (4) 

a'^hi{n) hi{n) > 1/a 



which is sketched in Fig. 0(a). We see that there is a hne of marginally stable fixed points 
hi G [a, 1). These fixed points are periodic states with period 1 — a: hi and /i2 take turns to 
topple and the toppling of one site will not trigger the toppling of another {hi < 1). Now, 
we introduce a little inhomogeneity. We drive hi with rate 1, but /i2 with a slightly slower 
rate (1 + e)~^. The Poincare map now reads 

hi{n) + e(l — a) a < hi{n) < 1 

hi{n + l) = \ 1 + a + e - a(l + t)hi{n) 1 < hi{n) < l/a (5) 

a^hiin) hiin) > l/a 

which is sketched in Fig. 0(b). In this case, there is only one fixed point /i* = 1 + yrf e- This 
fixed point is the phase-locked or synchronized state: the toppling of /i2 will trigger hi to 
topple {hi > 1). Note that Eq. (^) is only e away from Eq. (^), so the locking is rather weak 
and fiuctuations can play a crucial role. In the OFC model with open boundary conditions, 
the boundaries introduce inhomogeneity. The sites at the open boundaries have only three 
neighbors and hence have a slower effective growth rate. This inhomogeneity in the effective 
growth rate propagates into the interior of the sample, causing phase-locking and thus long 
range correlation. However, the whole system is not in a synchronized state. Rather, it is 
only "marginally" locked so that it gives a power-law distribution of avalanche sizes. 

We do not have a complete theory for the emergence of the "marginal locking" in the 
OFC model. However, we can abstract some of the features to construct a simpler model 
which also exhibits SOC This model is defined on a directed lattice to simplify avalanches. 
We define dynamical variables < 0j < 1 as the phase of the next toppling time of a site; this 
phase is related to the height hi in the OFC model at a fixed time, modulo the natural period 
(1 — 4a). At random sites on the boundary, we initiate a toppling which changes the phase 
according to (pi —>■ (p'^ = (f)i + a. If the phase (pj at a neighboring "downhill" site j = i + x or 
i + y (see inset in Fig. ^) meets the locking condition {(pi — (pj mod 1) < a, the site j locks 
onto the boundary site, with (pj — > 0^. This disturbance can then continue to propagate 
by sites further from the boundary becoming locked, in zero time. We refer to this locking 
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as "marginal" since neighboring phases only lock upon a crossing in the toppling time; 
the configuration is neutrally stable with respect to a continuous set of perturbations. This 
model differs from the OFC model most notably in the directed lattice and the instantaneous 
locking (even though avalanches in OFC occur in zero time, locking takes place over a time 
of order 1). 

A snapshot of the (pi in the steady state is depicted in Fig. §(a). The domains, bounded by 
solid lines, are regions where the toppling times are identical. Fig. 0(b) shows a configuration 
after an avalanche, with the dotted lines showing the previous domain configuration. An 
avalanche crosses domain boundaries only when the neighboring domains have times that 
differ by no more than a. Numerically, it is found that the number of domains ^(5*) of 
size S" at a fixed time has a power law tail that is independent of alpha, n{S) ~ S^^ , with 
a = 1.495 ± 0.005. Yet the avalanche distribution has an a dependent exponent, with the 
probability of an avalanche of size s behaving as -P(s) ~ s~'^^°'^ (Fig. |]). Note that domains 
and avalanches are closely related, with avalanches defining domains. The exponent a must 
satisfy the bound a < 3/2; the fact that our numerical result saturates this inequality, to 
within numerical error, suggests that the domains have linear roughness, with the average 
width of a domain proportional to the length of the domain perpendicular to the boundary. 

The relationship between r and a can be approximately explained. Assuming that 
the toppling times of neighboring domains are independent variables, the probability of an 
avalanche, which takes place in one domain, incorporating a given neighboring domain is just 
a. The incorporation of a neighboring domain increases the avalanche size to the scale of the 
neighboring domain, which is typically larger than the original domain. Given a distribution 
of domain sizes n{S) ~ S^'^ and a scale independent probability a' oc a that a larger domain 
than the initial one chosen will become part of the avalanche, followed by a probability a' 
that an even larger domain will become part of the avalanche, etc., leads to a distribution 
of avalanches P{s) ^ s^^*-"\ with T(a) = a + a'(l — a). This is in agreement with the 
linear fit of Fig. g with a' = (1.20 ± 0.04)a; the fitted value of r(0) = a = 1.505 ± 0.005 
agrees with the value determined by the domain distribution. The varying exponent of the 



avalanche distribution is then explained by the probability a' for the avalanche size to move 
out further into the tails of the a-independent domain size distribution. 

In this paper we have examined how SOC can arise in a model without a conservation 
law. The OFC model with periodic boundary condition has a continuous set of neutrally 
stable periodic states. In general, inhomogeneity destroys these periodic states and causes 
phase-locking which is the building block for long range correlations. We found that an open 
boundary results in the invasion of the interior by a marginally locked region in which the 
avalanche size distribution is a power-law. A simplified model on a directed lattice has been 
used to demonstrate how an a-dependent avalanche size distribution exponent can arise in 
such non-conserved dynamical models. Finally, we note that the OFC model is similar to 
the coupled "integrate-and-fire" oscillators studied in the context of neural networks and 
biology. A close cousin is the Peskin's model for the cardiac pacemaker. |T^ We found 



that the model in 2d with nearest neighbor coupling tends to lock into some periodic or 
"quasiperiodic" cluster state and that it has richer behaviors than simple synchronization. 
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FIGURES 
FIG. 1. Configurations of tlie Olami-Feder-Christensen model at various times after random 

initialization, demonstrating the invasion of the short-range correlated interior by the "SOC re- 
gion". The density of the color corresponds to the height variables < /ij < 1. The boundaries 
are periodic in the vertical direction and open in the horizontal. The lattice consists of 64? sites, 
with a = 0.07. Times are t = (a) 1.2 x 10^, (b) 2.4 x 10^, (c) 6.0 x 10^, (d) 36.0 x 10^. 

FIG. 2. A plot of the invasion distance vs. time for conservation parameters a = 0.07,0.15. 
Power law fits are shown as dashed lines; symbols indicate system size (A = 64^, □ = 128^, 
• = 256^). 

FIG. 3. Return map for the two-site system: (a) homogeneous system; (b) inhomogeneous 
system. 

FIG. 4. Domains in the simplified SOC model described in the text, for a 200^ directed 
lattice, (a) Configuration with domains of identical toppling times indicated by solid lines, (b) 
Configuration after an avalanche, with previous domains indicated by dashed lines. The inset 
sketches the directed lattice. 

FIG. 5. The avalanche distribution exponent T{a) for the directed model as a function of the 
conservation parameter a. The linear fit is discussed in the text. 
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